Introduction
Example of running model
- Running the model to equilibrium
- See the model behavior graphically
Add hypothetical disturbance
Effect of forest management on recovery resilience
Different management scenarios
- Relative variation: parameter variation using the original values
  - Different scenarios

Introduction

Here I am using a four state transition model (STM) to understand how forest management can increase forest resilience. The model has four process-parameters that define the transition rate between the four states: Regeneration, Boreal, Temperate and Mixed (Figure 1).

Each parameter represents a different process: $\alpha$ is the colonisation rate (R -> B, T, M), $\beta$ the succession (B, T -> M), $\theta$ the exclusion extintion (M -> B, T) and $\epsilon$ the pertubation (B, T, M -> R).

Example of running model

# call the model
source("vissault_model_v3.R")

Running the model to equilibrium

# parameters
params = read.table("pars.txt", row.names = 1)

# Wrapper to collect parameters for a given set of environmental conditions
pars = get_pars(ENV1 = 0, ENV2 = 0, params, int = 3)
pars

##      alphab      alphat       betab       betat       theta      thetat 
## 0.999991086 0.998278766 0.269252962 0.442427327 0.159791566 0.973086903 
##        epsB        epsT        epsM 
## 0.009397677 0.009397677 0.009397677

# Get the transition matrix
transition_matrix = get_matrix(ENV1 = 0, ENV2 = 0, int = 3)
transition_matrix$MAT

##              R            B           M           T
## R 9.906906e-01 0.0000598863 0.003873647 0.005375876
## B 3.885079e-05 0.9354824254 0.064478724 0.000000000
## M 3.896877e-03 0.0017832562 0.929843330 0.064476536
## T 5.373681e-03 0.0000000000 0.001804291 0.992822028

# Equilibrium of the model for the specific parameters (pars)
eq = get_eq(pars)
eq

## $eq
##           B           T           M 
## 0.004134084 0.571809507 0.414663890 
## 
## $ev
## [1] -0.1122421

$eq is the proportion of each state at equilibrium and $ev the eigenvalue.

See the model behavior graphically

Starting with 25% of the proportion for each state, we can graphically see the behavior of the model until reaching equilibrium.

# data frame
time <- seq(1, 50, by = 0.1)
dm2 <- data.frame(matrix(NA, nrow = length(time), ncol = 4))

# Initial condition
y = c(B = 0.25, T = 0.25, M = 0.25)
dm2[1, ] = c(y, 0.25)

# Function to get behavior
behavior <- function(ENV1, ENV2) {

  # paramters depending on environment
  pars = get_pars(ENV1 = ENV1, ENV2 = ENV2, params, int = 3)

  # getting eq
  for(i in 2:length(time)) {
    eq = runsteady(y = y, func = model, parms = pars, times = c(0, i))[[1]]
    eq[4] <- 1 - eq[1] - eq[2] - eq[3]
    dm2[i, ] <- eq
  }

  # plot
  par(family = 'serif', cex = 0.8, mar = c(4, 4, 1, 1))
  plot(time, dm2[, 1],
       type = "l",
       lwd = 1.5,
       ylim = c(0, 1),
       col = 2,
       ylab = "",
       xlab = "")
  lines(dm2[, 2], col = 3, lwd = 1.5)
  lines(dm2[, 3], col = 4, lwd = 1.5)
  lines(dm2[, 4], col = 5, lwd = 1.5)
}

Ev1 <- c(-1.75, -1.5, -1.25, -1, -0.75, -0.5, 0, 0.5, 1)
par(mfrow = c(3, 3))
for(i in 1: length(Ev1)) {
  behavior(ENV1 = Ev1[i], ENV2 = 0)
  mtext(paste("ENV1 = ", Ev1[i], sep = ""), 3, cex = 0.85)
  if(i == 4) mtext("State proportion", 2, line = 2, cex = 0.8)
  if(i == 1) legend(35, 0.45, c("B", "T", "M", "R"), col = 2:5, lty = 1, bty = "n", cex = 0.8)
}

Add hypothetical disturbance

Setting up the data frame with two different environmental conditions is a way to see the behavior of the model after a disturbance. In this case, the model will start with 25% of the proportion for each state and after reaching equilibrium, the environmnetal patters will change.

# Function to produce the data frame with the equilibrium proportion of each state
# based on the environmental variation.
behavior <- function(envComb1, envComb2) {
    xx <- seq(1, 40, by = 1)
    dat <- data.frame(matrix(NA, nrow = length(xx), ncol = 6))
    names(dat) <- c("ENV1", "ENV2", "B", "T", "M", "R")
    dat[c(1:20), c(1,2)] <- envComb1 #environmental 1
    dat[c(21: dim(dat)[1]), c(1,2)] <- envComb2 #environmental 2

    # Initial condition
    dat[1, c(3:6)] = c(B = 0.25, T = 0.25, M = 0.25, R = 0.25)

    # Behavior
    for(i in 2:length(xx)) {
        pars = get_pars(ENV1 = dat[i, 1],ENV2 = dat[i, 2], params, int = 3)
        y = c(B = dat[(i - 1), 3], T = dat[(i - 1), 4], M = dat[(i - 1), 5])
        eq = runsteady(y = y, func = model, parms = pars, times = c(0, i))[[1]]
      eq[4] <- 1 - eq[1] - eq[2] - eq[3]
      dat[i, c(3:6)] <- eq
    }
    return(dat)
}

# Run the function
dat <- behavior(envComb1 = c(-.05, -.05), envComb2= c(0.05, 0.05))

# plot
par(family = 'serif', cex = 0.8, mar = c(4, 4, 1, 1))
plot(0,
     xlim = c(0, dim(dat)[1]),
     ylim = c(0, max(dat[,c(3:6)]) + 0.1),
     xlab = "time",
     ylab = "state proportion"
     )
lines(dat$B, col = 2, lwd = 1.5)
lines(dat$T, col = 3, lwd = 1.5)
lines(dat$M, col = 4, lwd = 1.5)
lines(dat$R, col = 5, lwd = 1.5)
legend(34, 0.25, c("B", "T", "M", "R"), col = 2:5, bty = "n", lty = 1)

Effect of forest management on recovery resilience

Using a numerical approach, I will set forest management into the model by changing the parameters related to each management process. For example, plantation can enhance colonization rate ($\alpha$), pre-commercial thinning can enhance competitive exclusion ($\theta$), enrichment planting can enhance succession rate ($\beta$), and cutting can enhance disturbance rate ($\epsilon$; Figure 2). Dotted line is the base natural process that occur without intervention.

By simply increasing the value of each parameter related to management, we can simulate the inclusion of management in the model processes (see figure above). The recovery resilience, or time rate in which a system returns to equilibrium after a disturbance, is measured by the largest real part of the eigenvalue. The eigenvalue is optained by the Jacobian matrix and a nice example can be found in this vignette of the rootSolve package.

The choice of parameters is delicate in this kind of “sensitivity analysis”, then I chose to do 2 different tests in the parametric variation to simulate the response of resilience to the increasing in forest management.

Test1 varies a specific parameter (the tested one) from 0 to 1, and the other parameters remains with the original values (fitted with field data).

Test2 varies a specific parameter from 0 to 1, and the other parameters remains with a constant value constPar. So we tested with four constant values: 0.2, 0.4, 0.6, 0.8.

Test 1 - paramenters constant with original value

# Constant value for the main parameters from 0 to 1 + original value for the other parameters
# And environment 1 variation
int <- 3
parSeq <- seq(0.01, 1.1, 0.05)
Ev1 <- c(0, -1.5, -0.75, 0.5)

# running eigenvalue to each parameter
eql <- as.list("NA")
df <- data.frame()
  for(k in 1: length(pars)) {
    for(i in 1: length(Ev1)) {
        pars = get_pars(ENV1 = Ev1[i], ENV2 = 0, params, int = int)
        for(j in 1: length(parSeq)) {
            pars[k] = parSeq[j]
            df[j, 1]    <- parSeq[j]
            df[j, (i + 1)] <- get_eq(pars)$ev
        }}
eql[[k]] <- df
}

Simulations with ENV1 = 0 (Temperate for about Montreal latitude)

Pars <- c(expression(alpha), expression(alpha), expression(beta),
          expression(beta), expression(theta), expression(theta),
          expression(epsilon), expression(epsilon), expression(epsilon))
st <- c("B", "T", "B", "T", "", "T", "B", "T", "M")
par(family = 'serif', cex = 1.2, mfrow = c(3,3), mai = c(0.3, .5, .2, .2))
for(i in 1: length(pars)) {
    plot(eql[[i]]$V1, eql[[i]]$V2, type = "l", lwd = 1.5, xlab = "", ylab = "", ylim = c(-.45,0), xlim = c(0, 1.1))
  abline(v = get_pars(ENV1 = 0, ENV2 = 0, params, int = int), col = "gray55", lty = 3)
  if(i == 4) mtext(side = 2, "largest real part", line = 2.1, cex = 1.2)
  if(i == 2) mtext(Pars, at = get_pars(ENV1 = 0, ENV2 = 0, params, int = int), line = 0, cex = 0.75)
  legend("bottomleft", Pars[i], bty = "n")
  legend(0.005, -.41, st[i], bty = "n")
}

Simulations for different environments (-1.5 for Boreal, -0.75 for Mixed and 0.5 for Temperate)

par(family = 'serif', cex = 1.2, mfrow = c(3,3), mai = c(0.3, .5, .2, .2))
for(i in 1: length(pars)) {
    plot(eql[[i]]$V1, eql[[i]]$V3, type = "l", lwd = 1.5, xlab = "", ylab = "", ylim = c(-.5,0), xlim = c(0, 1.1))
  lines(eql[[i]]$V1, eql[[i]]$V4, lwd = 1.5, lty = 2)
  lines(eql[[i]]$V1, eql[[i]]$V5, lwd = 1.5, lty = 3)
  if(i == 4) mtext(side = 2, "largest real part", line = 2.1, cex = 1.2)
  if(i == 6) legend("bottomright", c("B (-1.5)", "M (-0.75)", "T (0.5)"), lty = 1:3, lwd = 1.5, bty = "n", cex = 1.2)
  legend("bottomleft", Pars[i], bty = "n")
  legend(0.01, -.45, st[i], bty = "n")
}

Interaction between parameters

Some practices can change more than a parameter, so here I want to see the interaction between these parameters influenced by the same practices. The follow function can generate a plot to visualize the interaction between $\alpha$ Boreal and $\alpha$ Temperate, $\beta$ Boreal and $\beta$ Temperate and $\theta$ and $\theta$ Temperate.

# function to generate the plot interaction between parameters
parInt <- function(Par1, Par2, ENV1) {

    # parameter variation
    int <- 3
    par1 <- seq(0.01, 0.99, 0.05)
    par2 <- seq(0.01, 0.99, 0.05)

    # all possible combination
    BT <- expand.grid(p1 = par1, p2 = par2)

    # get eigenvalue for each combination
    for(i in 1:dim(BT)[1]) {
      pars = get_pars(ENV1 = ENV1, ENV2 = 0, params, int = int)
      pars[c(Par1, Par2)] <- c(BT$p1[i], BT$p2[i])
      BT[i, 3] <- get_eq(pars)$ev
    }

    # setting eigenvalue dependence color
    rbPal <- colorRampPalette(c('red','blue'))
    BT$Col <- rbPal(20)[as.numeric(cut(BT$V3,breaks = 15))]

    # plot
    library(grDevices)
    par(family = "serif", cex = 0.8, mar = c(4, 4, 1, 1))
    layout(matrix(1: 2, ncol = 2), width = c(2, 1), height = c(1, 1))
    plot(BT$p1, BT$p2, pch = 15, cex = 2.9, col = BT$Col, xlab = names(pars[Par1]), ylab = names(pars[Par2]))
    mtext(paste("ENV1 = ", ENV1, sep = ""), 3, cex = 1)

    # legend
    legend_image <- as.raster(matrix(rbPal(12), ncol=1))
    plot(c(0, 2),c(0, min(BT$V3)), type = 'n', axes = F, xlab = '', ylab = '', main = 'eigenvalue')
    text(x = 1, y = seq(max(BT$V3), min(BT$V3), l = 5), labels=round(seq(min(BT$V3), max(BT$V3),l=5),3))
    rasterImage(legend_image, min(BT$V3), min(BT$V3), 0, 0)
}

$\alpha$B and $\alpha$T interaction

parInt(1, 2, ENV1 = -1.5)

parInt(1, 2, ENV1 = -.75)

parInt(1, 2, ENV1 = .5)

$\beta$B and $\beta$T interaction

parInt(3, 4, ENV1 = -.75)

$\theta$ and $\theta$T interaction

parInt(5, 6, ENV1 = -1.5)

parInt(5, 6, ENV1 = -.75)

parInt(5, 6, ENV1 = .5)

$\epsilon$M and $\epsilon$T interaction

parInt(9, 8, ENV1 = -1.5)

parInt(9, 8, ENV1 = -.75)

parInt(9, 8, ENV1 = .5)

Test 2 - parameters constant with the same value

# Constant value for all parameters from 0 to 1
int <- 3
parSeq <- seq(0.01, 0.99, 0.05)
constPar <- c(0.2, 0.4, 0.6, 0.8) #Constant value of all other parameters

# running eigenvalue to each parameter
pars = get_pars(ENV1 = 0, ENV2 = 0, params, int = int)
eql1 <- as.list("NA")
df <- data.frame()
for(k in 1: length(pars)) {
    pars = get_pars(ENV1 = 0, ENV2 = 0, params, int = int)
  for(l in 1: length(constPar)) {
    pars[-k] <- constPar[l]
    for(j in 1: length(parSeq)) {
        pars[k] = parSeq[j]
        df[j, 1]    <- parSeq[j]
        df[j, l + 1] <- get_eq(pars)$ev
     }
  }
eql1[[k]] <- df
}

# plot
source("/Users/wvieira/GitHub/NativeFunctions/addTransColor.R") #get transparance function
Pars <- c(expression(alpha), expression(alpha), expression(beta),
          expression(beta), expression(theta), expression(theta),
          expression(epsilon), expression(epsilon), expression(epsilon))
st <- c("B", "T", "B", "T", "", "T", "B", "T", "M")
par(family = 'serif', cex = 1.2, mfrow = c(3,3), mai = c(0.3, .5, .2, .2))
for(i in 1: length(pars)) {
    plot(0, type = "n", lwd = 1.5, xlab = "", ylab = "", xlim = c(0, 1), ylim = c(-0.35, 0.01))
  if(i == 4) mtext(side = 2, "largest real part", line = 2.1, cex = 1.2)
  if(i == 9) legend("bottomright", legend = constPar, lty = 1, col = c(1:4), bty = "n", cex = 1.2)
  legend("bottomleft", Pars[i], bty = "n")
  legend(0.01, -.32, st[i], bty = "n")
      points(eql1[[i]]$V1, eql1[[i]]$V2, type = "l", lwd = 1.5)
      points(eql1[[i]]$V1, eql1[[i]]$V3, type = "l", lwd = 1.5, col = 2)
      points(eql1[[i]]$V1, eql1[[i]]$V4, type = "l", lwd = 1.5, col = 3)
      points(eql1[[i]]$V1, eql1[[i]]$V5, type = "l", lwd = 1.5, col = 4)
    #add vertical lines
    for(j in 1:4) abline(v = constPar[j], col = addTrans(j, 120), lty = 3)
}

Testing all parameters combination

Now we saw the variation of each parameter, we want to see every parameter changing with each other to test all combinations of parameters. The sensitivity analysis helps us to test some of the following questions:

How is the relation between input and outputs of the model?
For what parameter are the outputs varying more?
Can we reduce uncertainty and hence increase robustness?

All these questions will then be guided to measure the impact of forest management on the parameters (based in figure 2).

Using the expand.grid function, the parameters variation goes from 0.09 to 0.99:

# parameters
parm <- seq(0.09, 0.99, 0.15)

# all possible combination
sim <- expand.grid(alphab = parm,
                   alphat = parm,
                   betab = parm,
                   betat = parm,
                   theta = parm,
                   thetat = parm,
                   eps = parm)

dim(sim)

## [1] 823543      7

# get eigenvalue for each variation (time ~ 29 hrs in one core):
# system.time(for(i in 1: dim(sim)[1]) {
#   eq <- get_eq(sim[i, ])
#   sim[i, 8] <- eq$ev
#   sim[i, 9] <- eq$eq[1] #eq B
#   sim[i, 10] <- eq$eq[2] #eq T
#   sim[i, 11] <- eq$eq[3] #eq M
#   cat(100*i/dim(sim)[1], "%", "\r")
# })

# save the output
# save(sim, file = "data/ev.RData")

# load the output
load("data/evq.RData")

Exploratory analysis

I’m going to start with an analysis of variance type III to get what parameter influenced more the eigenvalue.

library(car)

# Call anova eigenvalue (ev ~ all parameters variation + interactions)
aov <- Anova(lm1 <- lm(V8 ~
                        alphab + alphat + betab + betat + theta + thetat + eps +
                        alphab * alphat + alphab * eps + alphat * eps + betab * theta +
                        betab * thetat + betat * theta + betat * thetat,
                data = sim),
                singular.ok = FALSE,
                type = "III"
)

This is

# table
knitr::kable(aov)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	174.988458	1	10752.46026	0
alphab	114.598061	1	7041.67069	0
alphat	168.456032	1	10351.06435	0
betab	154.756488	1	9509.27281	0
betat	19.017266	1	1168.54788	0
theta	94.978393	1	5836.10716	0
thetat	1.290286	1	79.28377	0
eps	6239.745370	1	383411.65399	0
alphab:alphat	182.431943	1	11209.83773	0
alphab:eps	2567.238242	1	157748.27373	0
alphat:eps	2761.935300	1	169711.76205	0
betab:theta	75.175726	1	4619.29899	0
betab:thetat	40.836284	1	2509.25420	0
betat:theta	56.478228	1	3470.39978	0
betat:thetat	31.226060	1	1918.73779	0
Residuals	13402.318296	823528	NA	NA

This is

# function to calculate Omega²
source("/Users/wvieira/GitHub/NativeFunctions/CalcOmegaSq.R")
# get Omega²
om <- calcOmegaSq(aov, Order = TRUE)

# table
knitr::kable(om)

	Effect	omegaSq	omegaPercent
7	eps	0.2408186	49.8861617
10	alphat:eps	0.1065946	22.0813335
9	alphab:eps	0.0990804	20.5247453
8	alphab:alphat	0.0070402	1.4583995
2	alphat	0.0065008	1.3466631
3	betab	0.0059721	1.2371363
1	alphab	0.0044222	0.9160727
5	theta	0.0036650	0.7592149
11	betab:theta	0.0029007	0.6008941
13	betat:theta	0.0021791	0.4514090
12	betab:thetat	0.0015754	0.3263529
14	betat:thetat	0.0012045	0.2495199
4	betat	0.0007333	0.1519115
6	thetat	0.0000492	0.0101856

Different management scenarios

Here I will test different scenarios of management practices to understand their impact on forest resilience. Each scenario will be created following the same numerical approach explained in this section. Different from test 1, where we consider the increase in management practices for each parameter separated, here we will consider the interaction between practices because a management practice is never alone (e.g. after harvesting an area, plantation is expected; when thinning an area, future harvest and hence plantation is expected).

Considering the four main management practices (plantation, thinning, Enrichment plantation, and harvest), we will start with six different scenarios: (i) increasing harvest rate leads to increase in plantation rate; (ii) increasing thinning rate, leads to increase harvest and hence plantation rate. (iii) Increasing in enrichment leads to increase harvest and hence plantation, (iv) increasing in thinning and enrichment rate leads to increase harvest and hence plantation, (v) increasing thinning leads to increase harvest and (vi) just exploitation (harvest).

Scenarios:

plantation + harvest
plantation + thinning + harvest
plantation + enrichment + harvest
plantation + thinning + enrichment + harvest
thinning + harvest
harvest

Each of these two scenarios, will be multiplied by two for management focusing in the (i) Boreal forest or in the (ii) Temperate one.

These scenarios are summarized in the following table:

Parameters	Scenario 1	Scenario 1	Scenario 2	Scenario 2	Scenario 3	Scenario 3	Scenario 4	Scenario 4	Scenario 5	Scenario 5	Scenario 6	Scenario 6	Scenario 6
	B	T	B	T	B	T	B	T	B	T	B	M	T
$\alpha$B	↑	-	↑	-	↑	-	↑	-
$\alpha$T	-	↑	-	↑	-	↑	-	↑
$\beta$B	-				↑	-	↑	-
$\beta$T					-	↑	-	↑
$\theta$			↑	-			↑	-	↑	-
$\theta$T			-	↑			-	↑	-	↑
$\epsilon$B	↑	-	↑	-	↑	-	↑	-	↑	-	↑	-	-
$\epsilon$M	-	-	-	-	-	-	-	-	-	-	-	↑	-
$\epsilon$T	-	↑	-	↑	-	↑	-	↑	-	↑	-	-	↑

Using scenario 1 as an example, to increase harvest and plantation rate we will increase the parameters $\epsilon$ and $\alpha$, respectively. If the objective is to focus on Borel, we will increase $\alpha$B and keep $\alpha$T with the original (or natural) value.

Model outputs:

Other than evaluationg the largest real part of the eigenvalue as a resilience measure, I will evaluate two other outputs. First, the proportion of each state as a forest productivity meansure. And second, the transition probability (from the transition matrix) from T -> M and from M -> B. Using the transition probability I can estimate the migration rate towards North, where a bigger transition means bigger adaptation to climate change.

Relative variation: parameter variation using the original values

Varying the parameters relatively is a good way to take the original value of each parameter into account. However, this variation may not represent very well the effect of management practice (e.g. plantation can have a bigger impact than harvest because $\alpha$ is much bigger alpha = 0.9 than $\epsilon$ eps = 0.09. For that, I’m going to increase the impact of harvest on disturbance because the value of $\epsilon$ is too small comparing with other parameters.

# Sequence from 0 to 100% of Increase
sq <- seq(0.1, 1, 0.1) # 0 to 100% for all parameters except epsilon
sqH <- seq(0.3, 3, 0.3) # 30 to 300% for epsilon

# function to get scenarios
managSim2 <- function(Par1, ENV1, managLimit = NULL) {
  # data frame
  dat <- data.frame(x1 = get_pars(ENV1 = ENV1, ENV2 = 0, params, int = 3))
  for(i in 1: length(sq)) {
    dat[, i+1] <- dat$x1 + (dat$x1*sq[i])
  }

  # update disturbance (30 to 300%)
  for(i in 1: length(sqH)) {
    dat[c(7, 8, 9), i+1] <- dat$x1[7] + (dat$x1[7]*sqH[i])
  }

  # definitions
  if(!is.null(managLimit)) {
    manLim = managLimit
  }else manLim = length(sq)

  # data frame
  # X2 = Par1 eigenvalue; X3:X6 = Par1 state proportion
  egv <- data.frame(matrix(ncol = 6, nrow = manLim + 1))
  egv[1] <- c(0, (sq[1: manLim] * 100))
  pars <- get_pars(ENV1 = ENV1, ENV2 = 0, params, int = 3)

  # get eigenvalue for Par1
  for(i in 1: (manLim + 1)) {
    pars[Par1] <- dat[Par1, i]
    eq <- get_eq(pars)
    egv[i, 2] <- eq$ev #eigenvalue
    egv[i, c(3:5)] <- eq$eq # B, T, M
    egv[i, 6] <- 1 - sum(eq$eq) # R
  }

  # get transition probability
  if(ENV1 == -1.5) {
    a = 2; b = 3 # B -> M
    c = 2; d = 1; e = 1; f = 3 # B -> R -> M
    legd = c("B -> M", "B -> R -> M")
  }
  if(ENV1 == -0.75 | ENV1 == 0.5) {
    a = 3; b = 4 # M -> T
    c = 3; d = 1; e = 1; f = 4# M -> R -> T
    legd = c("M -> T", "M -> R -> T")
  }

  prob <- data.frame(matrix(ncol = 5, nrow = manLim +1))
  prob[1] <- c(0, (sq[1: manLim] * 100))
  pars <- get_pars(ENV1 = ENV1, ENV2 = 0, params, int = 3)
  eq <- get_eq(pars)[[1]]
  for(k in 1: (manLim + 1)) {
    pars[Par1] <- dat[Par1, k]
    eq <- get_eq(pars)[[1]]
    MAT <- get_matrix(ENV1 = ENV1, ENV2 = 0, pars = pars, eq = eq)$MAT
    prob[k, 2] <- MAT[a, b]
    prob[k, 3] <- MAT[c, d] * MAT[e, f]
  }

  # plot 1
  par(mfrow = c(1, 3), family = "serif", mar = c(4, 4, 1, 1), cex = 1.2)
  plot(egv$X1, egv$X2, type = "l", ylim = c(-0.27, -0.01), lwd = 1.75,
       xlab = "", ylab = "Largest real part")

  #plot 2
  plot(egv$X1, egv$X3, type = "l", ylim = c(0, 1), lwd = 1.75,
       xlab = "Foret management rate (%)", ylab = "State proportion")
  lines(egv$X1, egv$X4, lwd = 1.75, col = 2)
  lines(egv$X1, egv$X5, lwd = 1.75, col = 3)
  lines(egv$X1, egv$X6, lwd = 1.75, col = 4)

  legend("topleft", c("B", "T", "M", "R"), lty = 1, col = c(1:4), bty = "n", cex = 0.9)
  mtext(paste("ENV1 = ", ENV1, sep = ""), 3)

  #plot 3
  plot(prob$X1, prob$X2, type = "l", lwd = 1.75, ylim = c(0, 0.14),
       col = 3, ylab = "Transition probability", xlab = "")
  lines(prob$X1, prob$X3, lwd = 1.75, col = 2)

  legend("bottomleft", legd, lty = 1, col = c(3, 2), bty = "n", cex = 0.9)
}

Different scenarios

Considering the original value of parameters by increasing the parameters from 0 to 100%, we can see the real impact of forest management on recovery resilience.

Scenario 1: ↑ plantation and ↑ harvest

managSim2(Par1 = c(2, 8), ENV1 = -1.5)

managSim2(Par1 = c(2, 9), ENV1 = -0.75)

managSim2(Par1 = c(2, 9), ENV1 = 0.5)

Scenario 2: ↑ plantation, ↑ thinning (M) and ↑ harvest

managSim2(Par1 = c(2, 6, 8), ENV1 = -1.5)

managSim2(Par1 = c(2, 6, 9), ENV1 = -0.75, managLimit = 6) #limited to 6 because after that the model bugs

managSim2(Par1 = c(2, 6, 9), ENV1 = 0.5)

Scenario 3: ↑ plantation; ↑ enrichment (B and T) and ↑ harvest

managSim2(Par1 = c(2, 4, 8), ENV1 = -1.5)

managSim2(Par1 = c(2, 4, 9), ENV1 = -0.75)

managSim2(Par1 = c(2, 4, 9), ENV1 = 0.5)

Scenario 4: ↑ plantation; ↑ thinning (M); ↑ enrichment (B and T) and ↑ harvest

managSim2(Par1 = c(2, 6, 4, 8), ENV1 = -1.5)

managSim2(Par1 = c(2, 6, 4, 9), ENV1 = -0.75)

managSim2(Par1 = c(2, 6, 4, 9), ENV1 = 0.5)

Scenario 5: ↑ thinning (M); harvest

managSim2(Par1 = c(6, 8), ENV1 = -1.5)

managSim2(Par1 = c(6, 9), ENV1 = -0.75, managLimit = 6)

managSim2(Par1 = c(6, 9), ENV1 = 0.5)

Scenario 6: ↑ harvest

managSim2(Par1 = 8, ENV1 = -1.5)

managSim2(Par1 = 9, ENV1 = -0.75)

managSim2(Par1 = 9, ENV1 = 0.5)

Linking State Transition Models and forest management

Willian Vieira

June, 28 2017 - last update: August, 10 2017